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Abstract 

The intra-annual dynamics of wood formation, which involves the passage of newly produced cells through three suc- 
cessive differentiation phases (division, enlargement, and wall thickening) to reach the final functional mature state, 
has traditionally been described in conifers as three delayed bell-shaped curves followed by an S-shaped curve. Here 
the classical view represented by the 'Gompertz function (GF) approach' was challenged using two novel approaches 
based on parametric generalized linear models (GLMs) and 'data-driven' generalized additive models (GAMs). These 
three approaches (GFs, GLMs, and GAMs) were used to describe seasonal changes in cell numbers in each of the 
xylem differentiation phases and to calculate the timing of cell development in three conifer species [Picea abies 
(L.), Pinus sylvestris L., and Abies alba Mill.]. GAMs outperformed GFs and GLMs in describing intra-annual wood 
formation dynamics, showing two left-skewed bell-shaped curves for division and enlargement, and a right-skewed 
bimodal curve for thickening. Cell residence times progressively decreased through the season for enlargement, 
whilst increasing late but rapidly for thickening. These patterns match changes in cell anatomical features within a 
tree ring, which allows the separation of earlywood and latewood into two distinct cell populations. A novel statistical 
approach is presented which renews our understanding of xylogenesis, a dynamic biological process in which the 
rate of cell production interplays with cell residence times in each developmental phase to create complex seasonal 
patterns. 

Key words: Cambial activity, conifers, generalized linear and additive models (GLMs and GAMs), Gompertz functions (GFs), 
timing of cell development, tree ring, wood formation, xylogenesis. 



Introduction 

Despite the key role of wood in terrestrial biosphere and 
human societies, our understanding of the complex biologi- 
cal processes involved in its formation remains fragmentary 
(Plomion et al., 2001; Chaffey, 2002). In particular, the 
dynamics (timing, duration, and rate) of the wood forma- 
tion phases (cambial division, cell enlargement, and second- 
ary cell wall formation) remains poorly explored, despite 
this being of fundamental importance in determining 
the quantity (number of cells produced) and quality (cell 



anatomical characteristics) of the wood formed (Vaganov 
et al, 2011). 

In recent years, a growing number of studies have been 
devoted to wood formation monitoring, and these have pro- 
vided excellent data sets (see, for example, Rossi et al., 20066, 
2008). However, most of these studies have only focused on 
key aspects of the process, such as the onset and cessation of 
xylem formation (Rossi et al., 201 1), the maximal rate of cell 
production (Rossi et al., 20066), or more general intra-annual 
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growth patterns (Camarero et al, 2010; Rathgeber et al, 
2011; Cuny et al, 2012). Thus, a complete description of 
wood formation dynamics can only be found in isolated stud- 
ies (Wodzicki, 1971; Horacek et al, 1999). 

During wood formation, cells successively pass through 
three phases to reach a permanently mature state: (i) division; 
(ii) enlargement; and (iii) wall thickening (Wilson et al, 1966). 
The temporal succession of the differentiation phases under- 
gone by the cells forms a spatial succession of developmental 
zones along the tree ring. Drawing inspiration from hydro- 
logical system modelling, the three developmental zones and 
the mature zone can be viewed as four interconnected res- 
ervoirs, with cells 'flowing' from one reservoir to the next as 
they differentiate (Supplementary Fig. SI A available at JXB 
online). In conifers, the seasonal dynamics of the number 
of cells passing through the successive developmental zones 
has traditionally been described as three delayed bell-shaped 
curves followed by an S-shaped curve (Supplementary Fig. 
SIB) (Rossi et al, 2006a; Lupi et al, 2010; Cuny et al, 2012). 
This classical view is based on the assumption that the flow 
of the source (i.e. the rate of cell production) follows a bell- 
shaped curve through the season (Larson, 1994) and drives 
the dynamics of the subsequent reservoirs. 

Through the growing season, cells accumulate asymmet- 
rically, with an initial positive exponential phase followed 
by gradual slowing. Consequently, the Gompertz function 
(GF) — an asymmetric sigmoid growth function, where the 
rate is maximal at one -third of the asymptote before to slow 
down progressively — has been used to fit the data and calcu- 
late the timing of cell development (Supplementary Fig. SIC 
at JXB online) (Camarero et al, 1998; Rossi et al, 2003). This 
method allows variations in cell residence times in the differ- 
ent zones of xylogenesis to be calculated through the season, 
but only where these are close to linear (Rossi et al, 2006a; 
Deslauriers et al, 2003). 

Some pioneer studies have found that developing tracheids 
spend less time in enlargement and more time in thickening 
as the growing season proceeds (Whitmore and Zahner, 1966; 
Skene, 1969, 1972; Wodzicki, 1971; Kutscha et al, 1975). 
Such findings argue against the idea of repeated bell-shaped 
curves, which implies that cell residence times change in par- 
allel between the successive reservoirs through the season. 
Instead, they suggest that developing tracheids accumulate 
in the enlargement reservoir at the beginning of the season, 
when the increase in cell diameter is greatest, and then accu- 
mulate in the thickening reservoir at the end of the season, 
when the main focus is on building and lignifying their thick 
cell walls. 

This led to the proposal of the following hypothesis: the 
rate of cell transition from one development zone to the next 
and the cell residence times in these zones are not monotonous 
through a growing season, so that the resulting changes in the 
number of cells in the successive zones should be intrinsically 
more complex than the classical view; additionally, it should 
be possible to see the footprint of cell anatomical features in 
the intra-annual dynamics of the process that produced them, 
so that we can distinguish the transition between earlywood 
(large, thin-walled cells) and latewood (small, thick-walled 



cells) in the intra-annual pattern of enlarging and thickening 
cell numbers. 

To test this hypothesis, three statistical modelling 
approaches were applied to three conifer species (to ensure 
robustness) over 3 years of monitoring (to ensure general- 
ity): (i) the traditionally used parametric GF, which fits bell- 
and S-shaped distributions and represents the classical view; 
(ii) another parametric approach using generalized linear 
models (GLMs); and (iii) a semi-parametric approach using 
generalized additive models (GAMs), which are more flex- 
ible and are thus able to highlight structures in the data dis- 
tribution that differ from bell- or S-shaped curves, thereby 
allowing the classical view of wood formation dynamics to 
be challenged. 

Materials and methods 

Study site, xylem sampling, and sample preparation 

Dominant and healthy Scots pines (Pinus sylvestris L.), silver firs 
(Abies alba Mill.), and Norway spruces [Picea abies (L.) Karst.] (n-5 
per species) were selected in a mixed stand located in the Vosges 
Mountains (48°29'N, 7°09'E, and 643 m a.s.l.), in northeast France 
(Supplementary Table SI at JXB online). Microcores were collected 
weekly from the tree stems from April to November 2007-2009 and 
prepared in the laboratory to obtain transverse sections (5-10 um 
thick), which were stained with cresyl violet acetate and perma- 
nently mounted on glass slides. For detailed information about the 
study site, tree characteristics, sampling strategy, and sample prepa- 
ration, see Cuny et al. (2012). 

Microscopic observations of xylem cell differentiation 

Overall, 1450 anatomical sections were observed using an optical 
microscope (AxioImager.M2; Carl Zeiss SAS, France) under visible 
and polarized light at x 100^400 magnification to distinguish the cells 
in the different developmental zones. Cambial cells had a rectangu- 
lar shape, were surrounded by a thin primary wall, and had small 
radial diameters, whereas cells in the enlargement zone still had thin 
primary walls but had larger radial diameters. In contrast to cells 
in the cambial and enlargement zones, cells in the thickening zone 
developed a secondary wall that appeared bright under polarized 
light because of its particular arrangement of cellulose microfibrils 
(Abe et ah, 1997). Cresyl violet acetate staining, whereby cellulose 
stains purple and lignin stains blue (Kutscha et al, 1975), was used 
to follow the advancement of lignification, with cells in the thicken- 
ing phase exhibiting purple and blue walls, indicating that lignifica- 
tion was in progress, whereas mature cells had entirely lignified and 
thus had completely blue walls. 

Gompertz function, and generalized linear and additive models 

For each sample, the radial number of cells in the cambial (n c , see 
Supplementary Table S2 at JXB online for the list of variables used 
along with their notations and units), enlargement (« E ), thickening 
(n T ), and mature (n M ) zones was counted along three radial cell files 
(see Supplementary Table S3 for a description of the number of 
cells in each zone and their variation). The number of cells from the 
previous year was also counted on three radial files per sample and 
used to standardize the number of cells of the current year to reduce 
within-tree growth variability (Rossi et al, 2003). A dedicated func- 
tion of the package CAVIAR (Rathgeber, 2012) was used to apply 
this standardization using R statistical software (R Development 
Core Team, 2011). 

Following Rossi et al. (2003), GFs were fitted to the standardized 
number of mature cells, which is cumulative. 
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Thus, the GF is defined as: 



N(t)=A.e 



(1) 



where N(t) is the number of mature cells at time t; A is the upper 
asymptote parameter representing the final number of cells; |3 is the 
x-axis placement parameter that reflects the location of the origin; 
and k is the growth rate parameter that determines the spread of the 
curve along the time axis. 

However, in contrast to previous studies (Rossi et al., 2003), here 
it was considered a better practice to fit the first derivative of the 
Gompertz function (GF') directly to the standardized number of 
cells in the division, enlargement, and thickening zones. 

Thus, GF' is expressed as: 



N(t) = Ax.e 



(2) 



GF and GF' parameters were estimated using the nls function in R, 
which determines the non-linear (weighted) least-squares estimates 
of the parameters of a non-linear model (Bates and Chambers, 
1992). 

Additionally, two new approaches were used to model intra- 
annual wood formation dynamics based on GLMs and GAMs. The 
advance in regression analysis provided by GLMs and GAMs has 
been an important statistical development over the last 30 years. 
GAMs are semi-parametric extensions of GLMs, which are them- 
selves mathematical extensions of classical linear models (Wood, 
2006). 

The advantages of GLMs are that they define (i) a link function 
between a non-linear response variable and explanatory variables; 
and (ii) a probability distribution for non-normal errors (McCullagh 
and Nelder, 1983). This makes them particularly appropriate for 
modelling count data (in our case cell number variations) by using a 
log as a link function, and a Poisson distribution for errors. 

GAMs are GLMs in which the linear predictor depends, in 
part, on a sum of smooth functions of the predictors. GAMs 
are referred to as being data driven rather than model driven, 
because the data determine the nature of the relationship between 
the response variable and the set of explanatory variables rather 
than assuming a certain type of parametric relationship (Hastie 
and Tibshirani, 1986). The strength of GAMs lies in their ability 
to deal with highly non-linear and non-monotonic relationships 
between the response and the set of explanatory variables, which 
helps develop ecological models that better represent the underly- 
ing data, thereby increasing our understanding of biological sys- 
tems (Guisan et al, 2002). 

In this study, the weekly count of cells in each zone of wood for- 
mation was expressed as a function of the day of the year: 



difference between the model predicted values and the observations 
in the units of the response: 



MAE: 



(4) 



where y\ is the ith observed value, and y i is the corresponding 

model-fitted value. 

One problem with the MAE is that it depends on the scale of 
the response, meaning that it cannot be compared across series. 
Therefore, to compare the fittings between the different zones, in 
which the cell numbers differed, the mean absolute percentage error 
(MAPE) was also computed (Armstrong and Collopy, 1992): 



MAPE = £ 



(5) 



where y\ is the ith observed value, y i is the corresponding fitted 
value, and y is the mean of the n observed values. 

In addition, the modelling efficiency (EF), a relative measure 
of the goodness-of-fit that gives the proportion of variation in the 
response captured by the model, was calculated (Mayer and Butler, 
1993). This statistic is similar to the coefficient of determination (r 2 ), 
but can be used for non-linear models: 



2>-y) 2 



(6) 



EF is close to 1 when there is a good fit, and is <0 when the 
model has equal or lower predictive power than the mean of the 
observations. 

The Akaike information criterion (AIC), which penalized the 
goodness-of-fit of the model by the number of parameters, was also 
calculated (Akaike, 1973): 



AIC=-2.1og(L) + 2.A: 



(7) 



log(n)=ex+s(d)+e 



(3) 



where L refers to the maximized value of the likelihood function for 
the fitted model and k is the number of parameters. AIC is a relative 
measure: it tells nothing about how well a model fits the data in an 
absolute sense, but is used to compare models. In a set of models, 
the model that has the lowest AIC is selected because it provides 
the best compromise between the goodness-of-fit and the number 
of parameters. 



where n is the vector of the weekly count of cells in the considered 
zone, d is the vector of the corresponding day of the year, s is an 
unspecified smooth function, a is the intercept, and e is the error 
term. GAMs were fitted in R using the mgcv package (Wood, 2006). 

Gompertz functions (GFs and GFs'), GLMs, and GAMs were 
fitted for every year on each individual tree. The values of the fitted 
models were then averaged over the three studied years in order to 
calculate means representing the general wood formation dynamics 
of a species. 

Goodness-of-fit measurement 

The goodness-of-fit of the GFs, GLMs, and GAMs was assessed 
by computing the mean absolute error (MAE), as this is the most 
natural and most easily interpretable measure of model accuracy 
(Willmott and Matsuura, 2005). The MAE gives the mean absolute 



Residence duration 

Four numerical functions (S) were defined to represent the cumula- 
tive number of cells: 

Scetm (t)=n c (t)+n E (t)+n T (t)+n M (t) (8) 

Setm (t)=n E (t)+n T (t)+n M (t) (9) 

S TM (t)=n T (t)+n M (t) (10) 

S M (t)=n M (t) (11) 

where n represents functions (GFs, GLMs, or GAMs) that predict 
the number of cells in the cambial (n c ), enlargement (w E ), thickening 
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(n T ), and mature (« M ) zones at date t. Because these numerical func- 
tions are strictly increasing, their inverse function (S^ 1 ) can be used 
to compute cell timings. For a cell i, the dates of entrance into the 
enlargement (f E>i ), thickening (t Ti ), and mature (r Mi ) zones were 
computed as: 

t^W'ffl (12) 

ti^Snf 1 © (13) 

tM^V 1 ® (14) 

For each approach, the time spent by each cell i in the enlargement 
(d E ;) and thickening (d Ti ) zones, and the total duration of cell forma- 
tion (<f Fi ) were calculated as: 



dE,i — tfi— tgj 


(15) 


dTj-tMi _ t T i 


(16) 


dF,i = t M j-t E i 


(17) 



Transition rate 

For each approach, the rate of cell production by the cambial zone 
(r c t ) and the entry rate of cells into the enlargement (e E ,tX thickening 
(e Tjt ), and mature (e M t ) zones at day t were computed as: 

r C,t = S C ETM (t)-ScETM (t-1) (18) 
e E ,t=SETM(t)-S E TM(t-l) (19) 

e T , t =S TM (t)-S TM (t-l) (20) 
e M ,,=S M (t)-S M (t-1) (21) 

Results 

Comparative performance of models 

GAMs gave the best fits to the changes in cell numbers in 
the different zones of xylogenesis, with errors (MAE and 
MAPE) being, on average, twice as low as those computed 
for GLM and GF fittings (Table 1). Model efficiencies were 
also, on average, 30% higher for GAMs than for GLMs and 
GFs, showing that GAMs explained cell number variations 
better than the other models. Finally, AICs computed from 
GAMs were 30% lower than those computed from GF and 
GLM, indicating that GAMs offer the best trade-off between 
goodness-of-fit and parsimony. 



Changes in cell numbers in the zones of xylogenesis 

The number of cambial cells increased rapidly and culmi- 
nated early, at the beginning of June, for all three species 
(Fig. 1 A-C). This remarkable pattern was well reproduced by 
the GAMs, which showed bell-shaped curves skewed to the 
left (Fig. 2A-C). In contrast, GLMs and GFs were not able to 
preserve this pattern, delaying the occurrence of maximal cell 
numbers until later in the season and underestimating them 
(Table 2). Later in the season, the number of cambial cells 



gradually decreased and finally returned to approximately its 
initial value. Once again, GAMs reproduced this progressive 
decrease better than GFs and GLMs. 

In the enlargement zone, the number of cells also increased 
rapidly at the beginning of the growing season, in April. The 
number culminated at the very beginning of June, before 
decreasing progressively back to zero in September-October, 
indicating the end of tree growth (Fig. 1D-F). All three 
approaches were able to capture this pattern, showing bell- 
shaped curves skewed to the left (Fig. 2D-F), with curves 
obtained using GFs being close to those obtained using 
GAMs. GLMs, however, were unable to follow the rapid 
early increase in enlarging cells, again resulting in a shift of 
the maximum to later in the season (Table 2). 

In the thickening zone, the number of cells rapidly increased 
and reached an initial peak in the first half of June (Fig. 1G- 
I). The number then plateaued or slightly declined over 2-3 
weeks, and then increased again to reach a second peak later 
in the season, which was greater than the first and thus cor- 
responded to the maximum number of cells observed in this 
zone. In spruce, this maximum occurred at the beginning of 
August, which was 2 weeks earlier than in fir (mid-August) 
and 1 month earlier than in pine (beginning of September). 
The difference in magnitude between the two peaks differed 
between species, with the second late peak being 30% greater 
than the first peak in fir and spruce, but 100% greater in pine. 
Following this late maximum, the number of cells rapidly 
decreased until November, when no more cells remained in 
the thickening zone, indicating the end of biomass allocation 
to the stem. Only GAMs were able to capture the complex 
pattern of cell number changes in the thickening zone, clearly 
showing a bimodal curve skewed to the right (Fig. 2G-I). 
In contrast, the less flexible GFs and GLMs presented bell- 
shaped curves that exhibited their unique maximum at the 
place of the local minimum between the early and late peaks 
indicated by the raw data and GAMs (Figs 1G-I, 2G-I). 

From the beginning of June to November, there was a 
steady accumulation of cells in the mature zone (Fig. 1J-L). 
Only GAMs were able to detect two periods of faster accumu- 
lation separated by a period of slower increase (Fig. 2J-L). 

Cell development duration 

All three modelling approaches allowed calculation of the 
time a particular tracheid spends in each zone of xylogen- 
esis during the course of its differentiation (Supplementary 
Fig. S2 at JXB online). All methods agreed on the general 
trends in cell development duration through a tree ring, 
showing that cells spent less time in the enlargement zone 
(d E ) and more time in the thickening zone (d T ) from the 
beginning to the end of a ring (Fig. 3). Furthermore, since 
d E was always less than or equal to d T , the total duration of 
tracheid formation (c/ F ) increased from the beginning to the 
end of a ring. 

However, a more detailed look showed that the three 
approaches presented contrasting results. GFs and GLMs, 
which returned similar curves when fitting cell numbers in 
the different zones of xylogenesis, exhibited contrasting 
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Table 1. Mean absolute error (MAE), mean absolute percentage error (MAPE), model efficiency (EF), and Akaike information criterion 
(AIC) computed from the fittings of the Gompertz functions (GFs), and generalized linear and additive models (GLMs and GAMs) on 
the number of cells in the cambial (CZ), diameter enlargement (EZ), wall thickening (TZ), and mature (MZ) zones for Scots pine, Norway 
spruce, and silver fir (n=5 per species) during 2007, 2008, and 2009 (means). 



Zone 


Species 


MAE (cell) 




MAPE (%) 




EF (%) 






AIC 






GF 


GLM 


GAM 


GF 


GLM 


GAM 


GF 


GLM 


GAM 


GF 


GLM 


GAM 


CZ 


Pine 


0.75 


0.74 


0.39 


1 1 


1 1 


6 


43 


45 


83 


246 


248 


227 




Fir 


0.83 


0.82 


0.55 


g 


1 0 


7 


29 


30 


61 


247 


249 


232 




Spruce 


0.73 


0.71 


0.39 


11 


1 1 


6 


50 


52 


82 


245 


245 


228 


EZ 


Pine 


0.61 


0.68 


0.28 


32 


35 


14 


81 


73 


94 


226 


256 


160 




Fir 


0.81 


0.88 


0.31 


44 


46 


17 


68 


59 


93 


254 


333 


164 




Spruce 


0.41 


0.46 


0.17 


35 


38 


15 


83 


75 


95 


227 


198 


135 


TZ 


Pine 


1.72 


1.73 


0.71 


34 


34 


14 


67 


67 


92 


300 


519 


236 




Fir 


3.02 


3.08 


1.06 


30 


30 


11 


74 


74 


96 


338 


805 


292 




Spruce 


1.48 


1.4 


0.45 


36 


33 


10 


78 


78 


97 


293 


456 


198 


MZ 


Pine 


1.64 


1.25 


1.01 


14 


11 


9 


94 


95 


97 


302 


320 


238 




Fir 


2.85 


2.24 


1.44 


11 


9 


5 


97 


97 


99 


330 


442 


271 




Spruce 


1.47 


1.39 


0.84 


12 


11 


7 


96 


96 


98 


295 


312 


217 




U 100 150 200 250 300 100 150 200 250 300 100 150 200 250 300 




100 150 200 250 300 100 150 200 250 300 100 150 200 250 300 

Day of year Day of year Day of year 



Fig. 1. Number of cells in the cambial, enlargement, thickening, and mature zones of wood formation for Scots pine, silver fir, and 
Norway spruce. Lines represent the mean of five trees per species per year. The upper x-axis represents the months. 
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Pine Fir Spruce 




100 150 200 250 300 100 150 200 250 300 100 150 200 250 300 

Day of year Day of year Day of year 



Fig. 2. Gompertz functions (GFs), and generalized linear and additive models (GLMs and GAMs) fitted to the weekly cell count in the 
cambial, enlargement, thickening, and mature zones for Scots pine, silver fir, and Norway spruce (n=5 per species) over 3 years (2007, 
2008, and 2009). The upper x-axis represents the months. 

monotonous patterns, which were always convex for GFs durations. The lack of fit of GFs and GLMs resulted in a sys- 
and always concave for GLMs. In contrast, GAMs gave an tematic underestimation of cell development duration (30% 
accurate fit, capturing realistic patterns of cell development and 10%, respectively) compared with GAMs. 



Table 2. Maximum number (nj and date of occurrence of the maximum number (tj of cells in the cambial (CZ), diameter enlargement 
(EZ), and wall thickening (TZ) zones of wood formation, computed from the fittings of the Gompertz functions (GFs), and generalized 
linear and additive models (GLMs and GAMs) for Scots pine, Norway spruce, and silver fir (n=5 per species) during 2007, 2008, and 
2009 (means +SEs). 



Zone 


Species 


n % (cell) 






t x (day of year) 






GF 


GLM 


GAM 


GF 


GLM 


GAM 


CZ 


Pine 


7.6 + 0.3 


7.7 + 0.3 


8.5 + 0.3 


158 + 9 


156 + 8 


148 + 10 




Fir 


9.1 +0.3 


9.1+0.3 


10.0 + 0.5 


162 + 15 


161+15 


148 + 14 




Spruce 


7.8 + 0.2 


7.7 + 0.2 


8.8 + 0.3 


163 + 6 


162 + 4 


151 +3 


EZ 


Pine 


4.5 + 0.3 


4.4 + 0.3 


4.5 + 0.3 


153 + 2 


164 + 2 


151 +4 




Fir 


5.3 + 0.5 


5.1+0.4 


5.5 + 0.5 


150 + 4 


161 +3 


149 + 6 




Spruce 


4.8 + 0.3 


4.7 + 0.3 


4.6 + 0.3 


152 + 3 


162 + 3 


154 + 2 


TZ 


Pine 


8.7 + 0.9 


9.2 + 0.8 


10.7 + 1.2 


218 + 4 


218 + 3 


247 + 6 




Fir 


18.8 + 1.8 


20.6 + 2.1 


21.3 + 2.1 


208 + 5 


211+4 


229 ±9 




Spruce 


10.2 + 0.9 


11.6+1.0 


11.4 + 1.2 


198±4 


202 + 4 


216 + 7 
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Pine 
Tree ring % 

0 20 40 60 80 100 



Fir 

Tree ring % 



20 40 60 80 100 0 



Spruce 
Tree ring % 

20 40 60 80 100 




1 5 10 15 20 25 

Cell number 



20 40 60 

Cell number 



10 20 30 

Cell number 



Fig. 3. Cell residence durations in the enlargement (d E ) and thickening (d T ) zones, and total duration of formation (d F ) of tracheids along 
tree rings in Scots pine, silver fir, and Norway spruce as computed from Gompertz functions (GFs), and generalized linear and additive 
models (GLMs and GAMs). 



The GAMs showed that d E was longest in pine (13 ±4 
d), ~30% shorter in spruce (9 ± 3 d), and 50% shorter in fir 
(6 ±2 d). d E was also maximal for the first-formed cells of 
the year, ~20 d in pine, and 14 d in fir and spruce (Fig. 3). 
It decreased by approximately two-thirds along the tree 
ring (with final durations of 6 d in pine, 5 d in fir, and 4 d 
in spruce). However, while d E decreased progressively in pine 
and spruce, in fir it decreased rapidly in the first 30% of cells, 
and then remained stable around its minimal value for the 
following cells. The GAMs also showed that the ring could 
be separated into three distinct groups of cells according to 
d T (Fig. 3D-F). The first group corresponded to the first 40% 
of cells, for which d T remained at a steady minimum of ~20 d 
(regardless of species). The second group corresponded to the 
middle 40% of cells, which showed a strong increase in d T ; the 
first cells in this second group entered thickening at the end of 
June, when the number of cells in this zone started to increase 
again (Fig. 2). The third group corresponded to the last 20% 
of cells, which had the maximal d T (53 ± 1 d in pine, 51 ± 1 d 
in fir, and 45 ±2 d in spruce); the cells belonging to this third 
group entered thickening in August and/or September, when 



the number of cells in this zone was maximal. On average, 
d T was 10% shorter in fir and spruce than in pine (32 ± 14, 
3 1 ± 1 1 , and 36 ± 14 d, respectively). 

The longer d E and d T in pine cells meant that d ¥ was also 
20% longer in pine than in fir and spruce (49 ± 10, 38 ± 13, 
and 40 ±8 d, respectively). d ¥ followed a complex pattern 
resulting from the combination of d E and d T . In the first 40% 
of cells, d F decreased slightly (from 41 to 36 d in pine, 37 to 
30 d in spruce, and 35 to 24 d in fir). It then increased rapidly 
to reach its maximum at ~80% of a ring (61 d in pine, 57 d in 
fir, and 51 d in spruce). In the last 20% of cells, d F decreased 
slightly for a few days. 

On average, 20% of d E was explained by d E and 80% by 
d T , with little variation between species in these proportions 
(15% versus 85% in fir, 20% versus 80% in spruce, and 25% 
versus 75% in pine). However, these relative contributions 
varied greatly according to the position of the cells along the 
ring: the relative contribution of d E to d E in the first 40% of 
cells was ~50% in pine, 40% in spruce, and 30% in fir, but this 
fell to ~10% for the last 20% of cells in a ring, regardless of 
species. 



1990 | Cunyefa/. 



Entry rates into the different zones 

Cell entry rates into the successive zones of xylogenesis exhib- 
ited more complex seasonal patterns when calculated using 
GAMs than the classic bell-shaped curves depicted when 
using GFs or GLMs (Fig. 4). The cell division rates (and so 
also the entry rates) depended on the productivity of the spe- 
cies studied, with fir reaching ~0.7 cells d -1 , spruce reaching 
0.35 cells d" 1 , and pine remaining at ~0.25 cells d" 1 . 

Despite these contrasting levels, the three species presented 
similar changes in their rates of entry into the different zones 
of xylogenesis during the growing period. Moreover, the divi- 
sion rate, and the enlargement and thickening entry rates 
followed close to parallel patterns, only being staggered by a 
few days. These three rates presented a bimodal pattern: the 
first peak occurred at the end of May, around the time when 
the number of cells culminated in the enlargement zone and 
slightly before the first peak in the thickening zone; and the 
second peak occurred 40 d later, in the first half of July, when 
the number of cells in the enlargement zone was decreasing 
and at the beginning of the second peak in the thickening 
zone. In contrast, the entry rate into the mature zone followed 



a different pattern (Fig. 4J-L), reaching a clear maximum at 
the beginning of July, at the time of the small dip in the num- 
ber of cells in the thickening zone, and then decreasing but 
remaining fairly high, before reaching a second peak at the 
end of the season (October), at the same time as the strong 
final decrease in the number of cells in the thickening zone. 

Discussion 

Modelling wood formation dynamics using GFs, GLMs, 
and GAMs 

Numerous studies on the intra-annual dynamics of xylogen- 
esis have been conducted during the last decade, to gain a 
better understanding of the influence of climate change on 
wood formation (Gricar et al., 2011). The main method used 
to investigate xylogenesis involves repeated sampling of the 
stem through a growing season (Rossi et al, 2006a). However, 
growth is subject to large and sudden variations along and 
around the stem (Wodzicki and Zajaczkowski, 1970), which 
makes it difficult to decide whether differences between sam- 
ples represent the temporal changes in cambial activity (the 
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Fig. 4. Rate of cell production (r c ) by the cambial zone, and rates of entrance into the enlargement (e E ), thickening (e T ), and mature (e M ) 
zones of wood formation as computed from Gompertz functions (GFs), and generalized linear and additive models (GLMs and GAMs). 
The upper x-axis represents the months. 



biological signal of interest) or random within-tree particu- 
larities ('noise'). To enhance the desired biological signal, 
Wodzicki (1971) sampled a large number of trees (60 in 1967, 
216 in 1968) and obtained excellent results on the kinetics of 
tracheid development. However, such large sampling is excep- 
tional considering the amount of work and cost involved. 

In this study, three statistical approaches (GFs, GLMs, and 
GAMs) were tested to extract a meaningful biological signal 
from naturally noisy data sets (groups of five trees from three 
species over 3 years). It was found that GFs and GLMs accu- 
rately represented the general changes in the number of cells 
in the enlargement and mature zones, but were inaccurate 
concerning the cambial and thickening zones. Therefore, it is 
argued that GFs and GLMs could be used to compare gen- 
eral growth patterns between groups of trees — GFs can be 
used to extract meaningful biological parameters such as the 
final number of produced cells and the value and occurrence 
of the maximal rate of cell production (see, for example, 
Cuny et ah 2012) — but it is argued that they are unsuitable for 
accurately characterizing intra-annual dynamics of the wood 
formation process (zone entry rates and/or cell timings). 

In contrast, GAMs captured in detail the intrinsic com- 
plexity of intra-annual wood formation dynamics, making 
them the only sound approach that is currently available for 
accurately characterizing this process. Their use solves the 
basic but thorny issues of sampling date heterogeneity, sam- 
pling frequency variability, and missing data. GAMs allow 
trees and groups of trees (e.g. across species, sites, or years) 
to be compared or pooled. One possible problem, however, is 
that GAMs can produce results that are 'biologically incoher- 
ent' (e.g. negative zone entry rates), in which case the degree 
of smoothing must be increased until coherent results are 
obtained (Wood, 2006). 

Because of their high flexibility, GAMs could be adapted 
to a wide range of conditions where GFs are inappropriate, 
for example accounting for an additional resting period dur- 
ing the summer drought period in Mediterranean regions 
(De Luis et ah, 2007; Camarero et ah, 2010), or fitting more 
complex growth patterns as found in ring-porous deciduous 
trees (Michelot et ah, 2012). Moreover, mixed-effect general- 
ized additive models (GAMMs) could accurately account for 
the non-independence of weekly cell count data (Zuur et ah, 
2009), providing a more robust description of the system. 

Intra-annual wood formation dynamics 

Intra-annual wood formation dynamics is generally described 
as three delayed bell-shaped curves for the numbers of cells in 
the cambial (CZ), enlargement (EZ), and thickening (TZ) zones, 
followed by an S-shaped curve for the mature (MZ) zone (Rossi 
et ah, 2006a; Lupi et ah, 2010; Cuny et ah, 2012). However, this 
traditional view, which has been shaped by the use of GFs, did 
not stand the test of a GAM approach. Rather, it was demon- 
strated that, for conifers in temperate forests, the general pat- 
tern of intra-annual wood formation dynamics is composed of 
two slightly delayed left-skewed bell-shaped curves representing 
CZ and EZ, followed by a right-skewed bimodal curve repre- 
senting TZ, and a double sigmoid curve representing MZ. 
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The extensive data set of Wodzicki 's seminal work (1971) 
was used to plot the changes in intra-annual cell numbers 
in the enlargement, thickening, and mature zones for Scots 
pine trees grown in Poland (Supplementary Fig. S3 at JXB 
online). The patterns thus established were in full agreement 
with those highlighted in this study, with the curve for enlarg- 
ing cells being clearly shifted to the left and culminating early 
(mid-May), while the curve for thickening cells reached a pla- 
teau (June) and then increased again to reach a late maximum 
(September), and the accumulation of mature cells clearly 
followed a double sigmoid curve. In the 1960s, some authors 
also observed that the number of cells culminated early in the 
season for enlargement and late for thickening (Wodzicki and 
Peda, 1963; Whitmore and Zahner, 1966), which is supported 
by the present findings. 

The calculations showed that the enlargement duration 
(d E ) decreased from 2-3 weeks to <1 week as the season 
progressed, while thickening duration (d T ) increased from 
3 weeks to nearly 2 months. These values match those of 
Wodzicki (1971), whereas Horacek et ah (1999) found longer 
d E but similar d T . Both Wodzicki (1971) and Horacek et ah 
(1999) also described the same characteristic changes in d T : 
steady minimal durations for the first 40% of cells in a ring, 
a strong increase for the following 40%, and a new plateau 
around the maximum duration for the last 20%. Conversely, 
studies based on GFs have found pseudo-linear variations 
in the differentiation durations between the successive cells 
of a ring (Deslauriers et ah, 2003; Rossi et ah, 2006a) and, 
moreover, have calculated shorter d T , from 1-3 weeks for the 
first tracheids to 4-5 weeks for the last tracheids. It should be 
noted that these studies were performed in colder environ- 
ments (boreal or subalpine zone), which may explain the dis- 
crepancies, but pseudo-linear variations and underestimated 
values of d T were also obtained by using GFs. 

The developing xylem is a dynamic biological system 

The present results show that the developing xylem is a 
dynamic biological system, justifying the analogy made with 
hydrological system modelling. It has been demonstrated 
that the cell flow through the different zones of xylogenesis 
is controlled not only by the production rate of the cambium 
(the source), but also by the cell residence times in each 'res- 
ervoir' (Supplementary Fig. S4A at JXB online), creating 
the complex pattern described above (Supplementary Fig. 
S4B), and resulting in strong differences in the computa- 
tion of the timing of cell development (Supplementary Fig. 
S4C). The early culmination of cells in the enlargement res- 
ervoir can be related to both (i) the high flow from the source 
(demonstrated by the high division rate in the cambium and 
the high entry rates into the enlarging reservoir); and (ii) the 
long cell residence times in this reservoir at the beginning 
of the season. The enlargement reservoir then progressively 
emptied, because of the decreasing residence times in this 
reservoir. 

The early peak in cells in the thickening reservoir can also 
be related to the intense flow from the source at the beginning 
of the growing season: the many cells born in the cambium 
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passed quite rapidly through the enlarging reservoir, before 
entering the thickening reservoir, where they accumulated 
because of structurally longer residence times. After this first 
peak, the number of cells in the thickening reservoir remained 
quasi-stable for a while, indicating that the entry and exit rates 
were at equilibrium. At the end of summer, despite a rapid 
decrease in the entry rate, the thickening reservoir reached its 
maximum cell number, because of an abrupt increase in the 
cell residence times. 

The complex dynamics of the thickening reservoir was 
demonstrated by the accumulation of mature cells, which had 
two periods of rapid increase separated by a slower transition 
period. The first rapid accumulation echoed the intense cell 
flow generated by the cambium in the first part of the sea- 
son, which cascaded through the successive reservoirs. In sum- 
mer, the period of slower accumulation corresponded to the 
increasing residence times in the thickening reservoir, result- 
ing in a diminution of mature reservoir entry rate. In autumn, 
the last period of rapid accumulation of mature cells occurred 
during the final rapid emptying of the thickening reservoir 

Relationships between the intra-annual dynamics of 
xylogenesis and wood anatomy 

Conifer tree rings can be divided into two parts: earlywood, 
containing large-diameter, thin-walled tracheids produced 
at the beginning of the growing season; and latewood, con- 
taining narrow-diameter, thick-walled tracheids produced 
towards the end of the season. In the present data set, the 
first latewood cell — identified according to Mork's criterion 
(Denne, 1988) — was estimated to enter the thickening zone 
around the beginning of the second peak in thickening cells. 
Therefore, the bimodal curve describing the number of cells 
in the thickening zone can be interpreted as two juxtaposed 
bell-shaped curves, corresponding to two distinct populations 
of tracheids (earlywood and latewood) exhibiting contrast- 
ing behaviours. The higher proportion of latewood in pine 
(45%) compared with spruce (30%) and fir (27%) explains the 
greater difference in amplitude between the two peaks in the 
pine thickening zone. 

The observed patterns for cell differentiation durations are 
consistent with the anatomical profiles of conifer tree rings 
(Makinen et al., 2003; Park and Spiecker, 2005; Rathgeber 
et al, 2006). Indeed, the final diameter of a tracheid, which is 
acquired during the enlargement phase of its differentiation, 
progressively decreases along a ring, mimicking the changes in 
d E ; and the final wall thickness, which is acquired during the 
thickening phase, remains around a steady minimum in the 
first part of a ring, increases in the middle part, and reaches a 
steady maximum in the last part, mimicking the changes in d T . 
Moreover, the sizes of tracheids are proportional to the time 
they spent in the different zones of differentiation: in earlywood 
tracheids, diameters are two to three times larger and d E two 
to three times longer than in latewood tracheids; whereas in 
latewood tracheids, walls are two to three times thicker and d T 
two to three times longer than in earlywood tracheids. It is also 
striking that the changes in d F along a ring strongly resembled a 
classic intra-ring wood density profile (Rathgeber et al., 2006). 



Conclusion 

The present study shows that GAMs provide an improved 
approach for modelling intra-annual wood formation dynam- 
ics, which is generally represented by changes in the number 
of cells in the different development stages of xylogenesis: 
cambial (CZ), enlargement (EZ), wall thickening (TZ), and 
mature (MZ) zones. Indeed, only GAMs captured the com- 
plexity of intra-annual wood formation dynamics, making 
them suitable for further characterizing this dynamics by 
computing both the cell entry rates and the cell residence 
times in the different zones. 

It is argued that the developing xylem must be seen as a 
dynamic biological system in which the cell production rate 
interplays with the cell residence time in each zone, resulting 
in complex intra-annual patterns: two left-skewed bell-shaped 
curves for CZ and EZ, followed by a right-skewed bimodal 
curve for TZ and a double sigmoid curve for MZ. These 
patterns have the advantage of truly reflecting the data and 
explaining some anatomical features of a tree ring, separat- 
ing, for example, earlywood and latewood into two distinct 
cell populations. 

It is believed that this approach has great potential and 
will help to better understand wood formation. It allows a 
fine assessment of the dynamics of the processes that can be 
compared directly with observations; for example, a compari- 
son of cell residence times in each of the zones with tracheid 
anatomical characteristics shows that the change in cell size 
along a ring is closely related to the change in d E through the 
growing season, whereas the change in cell wall thickness is 
closely related to the change in d T . 

Supplementarydata 

Supplementary data are available at JXB online. 

Figure SI. Classical view of the wood formation dynamics. 

Figure S2. Sums of the cell numbers in the wood formation 
zones as simulated by the models. 

Figure S3. Data used by Wodzicki (1971). 

Figure S4. Updated view of the wood formation dynamics. 

Table SI. Characteristics of the monitored trees. 

Table S2. List of the variables. 

Table S3. Cell numbers and their variations in the zones of 
xylogenesis. 
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